home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / sph_poly.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  2.4 KB  |  110 lines
  1. /*
  2.  * ANSI C code from the article
  3.  * "Computing the Area of a Spherical Polygon"
  4.  * by Robert D. Miller
  5.  * in "Graphics Gems IV", Academic Press, 1994
  6.  */
  7.  
  8. /*
  9.  * to test, compile with "cc -DMAIN -o spherical_poly spherical_poly.c -lm"
  10.  * to create subroutine, compile with "cc -c spherical_poly.c"
  11.  */
  12.  
  13. #include <math.h>
  14. static const double
  15.     HalfPi= 1.5707963267948966192313,
  16.     Degree= 57.295779513082320876798;    /* degrees per radian */
  17.  
  18. double Hav(double X)
  19. /*  Haversine function: hav(x)= (1-cos(x))/2.  */
  20. {
  21.     return (1.0 -cos(X))/2.0;
  22. }
  23.  
  24. double SphericalPolyArea(double *Lat, double *Lon, int N)
  25. /*  Returns the area of a spherical polygon in spherical degrees,
  26.      given the latitudes and longitudes in Lat and Lon, respectively.
  27.      The N data points have indices which range from 0 to N-1. */
  28. {
  29. int  J, K;
  30. double
  31.     Lam1, Lam2, Beta1, Beta2,
  32.     CosB1, CosB2, HavA,
  33.     T, A, B, C, S, Sum, Excess;
  34.  
  35.     Sum= 0;
  36.     for (J= 0; J <= N; J++)
  37.     {
  38.     K= J+1;
  39.     if (J == 0)
  40.     {
  41.         Lam1= Lon[J];    Beta1= Lat[J];
  42.         Lam2= Lon[J+1];    Beta2= Lat[J+1];
  43.         CosB1= cos(Beta1);    CosB2= cos(Beta2);
  44.     }
  45.     else
  46.     {
  47.         K= (J+1) % (N+1);
  48.         Lam1= Lam2;        Beta1= Beta2;
  49.         Lam2= Lon[K];    Beta2= Lat[K];
  50.         CosB1= CosB2;    CosB2= cos(Beta2);
  51.     }
  52.  
  53.     if (Lam1 != Lam2)
  54.     {
  55.         HavA= Hav(Beta2-Beta1) +CosB1*CosB2*Hav(Lam2-Lam1);
  56.         A= 2*asin(sqrt(HavA));
  57.         B= HalfPi -Beta2;
  58.         C= HalfPi -Beta1;
  59.         S= 0.5*(A+B+C);
  60.         T= tan(S/2) * tan((S-A)/2) * tan((S-B)/2) * tan((S-C)/2);
  61.  
  62.         Excess= fabs(4*atan(sqrt(fabs(T))))*Degree;
  63.         if (Lam2 < Lam1) Excess= -Excess;
  64.  
  65.         Sum= Sum + Excess;
  66.     }
  67.     }
  68.     return fabs(Sum);
  69. }   /*    SphericalPolyArea. */
  70.  
  71. #ifdef MAIN
  72.  
  73. #include <stdlib.h>
  74. #include <stdio.h>
  75.  
  76. double SphericalPolyArea(double *Lat, double *Lon, int N);
  77.  
  78. double
  79.     SqMi= 273218.4, /* Square mi per spherical degree. */
  80.     SqKm= 707632.4; /* Square km per spherical degree. */
  81.  
  82. main()         /*     Spherical Polygon Area test routine */
  83. {
  84. double    Lat[100], Lon[100];
  85. double    P, Q, Area;
  86. int    K, LastPt;
  87.  
  88.     LastPt= -1;
  89.     printf("Enter Longitude Latitude. End with: 0 0.\n");
  90.     while (1)
  91.     {
  92.     scanf("%lf %lf",&P,&Q);
  93.     if ((P == 0.0) && (Q == 0.0))
  94.     break;
  95.     Lon[++LastPt]= P;  Lat[LastPt]= Q;
  96.     }
  97.  
  98.     for (K= 0; K <= LastPt; K++)
  99.     {
  100.     Lat[K]= Lat[K]/Degree;    Lon[K]= Lon[K]/Degree;
  101.     }
  102.  
  103.     Area= SphericalPolyArea(Lat, Lon, LastPt);
  104.     printf("     Area=%12.4lf sq mi., %12.4lf spherical degrees.\n",
  105.        Area*SqMi, Area);
  106.     return 0;
  107. }
  108.  
  109. #endif
  110.